% !Mode:: "TeX:UTF-8"
\section{快速傅里叶变换}

\subsection{问题}

\begin{enumerate}
 \item 总结书本上快速傅里叶变换内容，开发计算傅里叶变换程序。要求能显示幅度谱和相位谱
 \item 对于2-D傅里叶变换需要能做到0点居中
 \item 能正确地进行反变换计算
\end{enumerate}


\subsection[背景原理]{背景原理\footnote{因为仅涉及到数字图像处理问题，所以本节只分析离散傅里叶变换。同时由于篇幅问题，本节中的公式、定理、性质均不加证明的给出}}
\subsubsection{1-D离散傅里叶变换}


\begin{cdefi}
对一个连续函数$f(x)$等间隔采样可得到一个离散序列，设共采样了N个样，则这个离散序列可以表示为$\{f(0),f(1),f(2),\ldots,f(N-1)\}$。借助这种表达，并令$x$为离散实变量， $u$为离散频率变量，可将离散傅里叶变换对定义为：
\begin{align}
\label{e:dft}
&\mathbf{F}\{f(x)\}=F(u)=\sum_{x=0}^{N-1}f(x)e^{-j2\pi ux/N} \quad u=0,1,\ldots,N-1\\
&\mathbf{F}^{-1}\{F(u)\}=f(x)=\frac{1}{N}\sum_{u=0}^{N-1}F(u)e^{j2\pi ux/N} \quad x=0,1,\ldots,N-1
\end{align}
$F(u)$是复函数，可以写成：
\begin{equation}
F(u)=R(u)+jI(u) \quad \mbox{$R(u)$和$I(u)$分别为$F(u)$的实部和虚部。}
\end{equation}
\end{cdefi}
其中
\begin{align}
&|F(u)|=\sqrt{R^2(u)+I^2(u)}\\
&\phi(u)=\arctan[I(u)/R(u)]\\
&P(u)=|F(u)|^2=R^2(u)+I^2(u)
\end{align}
三式中幅度函数$|F(u)|$也称为$f(x)$的傅里叶{\bf 频谱},$\phi(u)$称为{\bf 相位角}， $P(u)$称为$f(x)$的{\bf 功率谱}。

\subsubsection{1-D快速离散傅里叶变换}

先将式($\ref{e:dft}$)写成：

\begin{equation}
\label{e:fft1}
F(u)=\sum_{x=0}^{N-1}f(x)W^{ux}_N
\end{equation}

其中：

\begin{equation}
\label{e:fftWN}
W_N=e^{\frac{-j2\pi}{N}}
\end{equation}

设N为2的正整数次幂，即：

\begin{equation}
N=2^n
\end{equation}

如令M为正整数，且：

\begin{equation}
\label{e:fftN2M}
N=2M
\end{equation}

将式($\ref{e:fftN2M}$)代入式($\ref{e:fft1}$) 得到：
\begin{equation}
F(u)=\sum_{x=0}^{2M-1}f(x)W^{ux}_{2M}=
\frac{1}{2}
\left[
\sum_{x=0}^{M-1}f(2x)W^{u(2x)}_{2M}+\sum_{x=0}^{M-1}f(2x+1)W^{u(2x+1)}_{2M}
\right]
\end{equation}

由式($\ref{e:fftWN}$)可知$W^{2ux}_{2M}=W^{ux}_{M}$,所以可将上式写为：

\begin{equation}
\label{e:fft2}
F(u)=
\frac{1}{2}
\left[
\sum_{x=0}^{M-1}f(2x)W^{ux}_{M}+\sum_{x=0}^{M-1}f(2x+1)W^{ux}_{M}W^u_{2M}
\right]
\end{equation}

定义：

\begin{align}
\label{e:fftEven}
F_{even}(u)&=\sum_{x=0}^{M-1}f(2x)W^{ux}_N \quad u=0,1,\ldots,M-1\\
\label{e:fftOdd}
F_{odd}(u)&=\sum_{x=0}^{M-1}f(2x+1)W^{ux}_N \quad u=0,1,\ldots,M-1
\end{align}

就可以将($\ref{e:fft2}$)简化为：

\begin{equation}
F(u)=\frac{1}{2}\left[  F_{even}(u)+ F_{odd}(u)W^u_{2M}\right]
\end{equation}

同理，由$W^{u+M}_{M}=W^u_M$和$W^{u+M}_{2M}=-W^u_{2M}$可得：

\begin{equation}
F(u+M)=\frac{1}{2}\left[ F_{even}(u)- F_{odd}(u)W^u_{2M} \right]
\end{equation}

表明一个N点的变换可以通过将原始表达式分成两半来计算。可以证明：通过上述算法，N点变换的复杂度为$N\log_2N$。具体16点序列的FFT计算步骤如下图\footnote{来自\url{http://www.texample.net/tikz/examples/radix2fft/}}所示：


%figure FFT
\tikzstyle{n}= [circle, fill, minimum size=4pt,inner sep=0pt, outer sep=0pt]
\tikzstyle{mul} = [circle,draw,inner sep=-1pt]
% Define two helper counters
\newcounter{x}\newcounter{y}
\begin{tikzpicture}[yscale=0.5, xscale=1.2, node distance=0.3cm, auto]
    % The strategy is to create nodes with names: N-column-row
    % Input nodes are named N-0-0 ... N-0-15
    % Output nodes are named N-10-0 ... N-10-15
    % Draw inputs
    \foreach \y in {0,...,15}
        \node[n, pin={[pin edge={latex'-,black}]left:$x(\y)$}]
              (N-0-\y) at (0,-\y) {};
    % Draw outputs
    \foreach \y / \idx in {0/0,1/8,2/4,3/12,4/2,5/10,6,7/14,
                           8/1,9,10/5,11/13,12/3,13/11,14/7,15}
        \node[n, pin={[pin edge={-latex',black}]right:$X(\idx)$}]
              (N-10-\y) at (7,-\y) {};
   % draw connector nodes
    \foreach \y in {0,...,15}
        \foreach \x / \c in {1/1,2/3,3/4,4/6,5/7,6/9}
            \node[n, name=N-\c-\y] at (\x,-\y) {};
    % draw x nodes
    \foreach \y in {0,...,15}
        \foreach \x / \c  in {1/2,4/5,7/8}
            \node[mul, right of=N-\x-\y] (N-\c-\y) {${\times}$};
    % horizontal connections
    % Note the use of simple counter arithmetics to get correct
    % indexes.
    \foreach \y in {0,...,15}
        \foreach \x in {0,1,3,4,6,7,9}
        {
            \setcounter{x}{\x}\stepcounter{x}
            \path (N-\x-\y) edge[-] (N-\arabic{x}-\y);
       }
    % Draw the W_16 coefficients
    \setcounter{y}{0}
    \foreach \i / \j in {0/0,1/0,2/0,3/0,4/0,5/0,6/0,7/0,
                            0/1,1/1,2/1,3/1,4/1,5/1,6/1,7/1}
    {
        \path (N-2-\arabic{y}) edge[-] node {\tiny $W^{\i\cdot\j}_{16}$}
                (N-3-\arabic{y});
        \stepcounter{y}
    }
    % Draw the W_8 coefficients
    \setcounter{y}{0}
    \foreach \i / \j in {0/0,1/0,2/0,3/0,0/1,1/1,2/1,3/1,
                         0/0,1/0,2/0,3/0,0/1,1/1,2/1,3/1}
    {
        \path (N-5-\arabic{y}) edge[-] node {\tiny $W^{\i\cdot\j}_{8}$}
              (N-6-\arabic{y});
        \addtocounter{y}{1}
    }
    % Draw the W_4 coefficients
    \setcounter{y}{0}
    \foreach \i / \j in {0/0,1/0,0/1,1/1,0/0,1/0,0/1,1/1,
                            0/0,1/0,0/1,1/1,0/0,1/0,0/1,1/1}
    {
        \path (N-8-\arabic{y}) edge[-] node {\tiny $W^{\i\cdot\j}_{4}$}
              (N-9-\arabic{y});
        \stepcounter{y}
    }
    % Connect nodes
    \foreach \sourcey / \desty in {0/8,1/9,2/10,3/11,
                                   4/12,5/13,6/14,7/15,
                                   8/0,9/1,10/2,11/3,
                                   12/4,13/5,14/6,15/7}
       \path (N-0-\sourcey.east) edge[-] (N-1-\desty.west);
    \foreach \sourcey / \desty in {0/4,1/5,2/6,3/7,
                                   4/0,5/1,6/2,7/3,
                                   8/12,9/13,10/14,11/15,
                                   12/8,13/9,14/10,15/11}
        \path (N-3-\sourcey.east) edge[-] (N-4-\desty.west);
    \foreach \sourcey / \desty in {0/2,1/3,2/0,3/1,
                                   4/6,5/7,6/4,7/5,
                                   8/10,9/11,10/8,11/9,
                                   12/14,13/15,14/12,15/13}
        \path (N-6-\sourcey.east) edge[-] (N-7-\desty.west);
    \foreach \sourcey / \desty in {0/1,1/0,2/3,3/2,
                                   4/5,5/4,6/7,7/6,
                                   8/9,9/8,10/11,11/10,
                                   12/13,13/12,14/15,15/14}
        \path (N-9-\sourcey.east) edge[-] (N-10-\desty.west);
\end{tikzpicture}


可以看到,再此过程中输出序列需要重新排序。这种重新排序有一定的规律，称之为{\bf 码位倒置}\footnote{《信号与线性系统》阎鸿森等编\;西安交通大学出版社\; 1999}。即先将序列的序号用二进制表示，然后将二进制码位倒置。再将倒置后的二进制码还原为十进制数，在该处放置相应序号的序列值。以N=16为例，码位倒置过程如下表：

\begin{longtable}{|c|c|c|c|c|c|c|c|c|}
\caption{码位倒置} \\
%\toprule
%\multicolumn{9}{c}{码位倒置} \\
%\midrule
\endfirsthead
%\midrule
%\multicolumn{9}{l}{码位倒置} \\
\midrule
\endhead
\midrule
\multicolumn{9}{l}{Continued\dots} \\
\endfoot
\endlastfoot
\bottomrule
序列值  & $x(0)$ & $x(1)$ & $x(2)$ & $x(3)$ & $x(4)$ & $x(5)$ & $x(6)$ & $x(7)$\\
\hline
序码号  & 0000   & 0001   & 0010   & 0011   & 0100   & 0101   & 0110   & 0111 \\
\hline
倒置    & 0000   & 1000   & 0100   & 1100   & 0010   & 1010   & 0110   & 1110 \\
\hline
重排序列& $x(0)$ & $x(8)$ & $x(4)$ & $x(12)$& $x(2)$ & $x(10)$& $x(6)$ & $x(14)$\\
\bottomrule
序列值  & $x(8)$ & $x(9)$ & $x(10)$& $x(11)$& $x(12)$& $x(13)$& $x(14)$& $x(15)$\\
\hline
序码号  & 1000   & 1001   & 1010   & 1011   & 1100   & 1101   & 1110   & 1111 \\
\hline
倒置    & 0001   & 1001   & 0101   & 1101   & 0011   & 1011   & 0111   & 1111 \\
\hline
重排序列& $x(1)$ & $x(9)$ & $x(5)$ & $x(13)$& $x(3)$ & $x(11)$& $x(7)$ & $x(15)$\\
\bottomrule
\end{longtable}

\subsubsection{2-D快速离散傅里叶变换}

\begin{cdefi}
将1-D离散傅里叶变换推广到2-D，有以下变换对
\begin{align}
\label{e:dft2d}
F(u,v)=\sum_{x=0}^{N-1}\sum_{y=0}^{N-1}f(x,y)e^{-j2\pi(ux+vy)/N}  \quad &u,v=0,1,\ldots,N-1\\
f(u,v)=\frac{1}{N^2}\sum_{u=0}^{N-1}\sum_{v=0}^{N-1}F(u,v)e^{j2\pi(ux+vy)/N}   \quad &x,y=0,1,\ldots,N-1
\end{align}
\end{cdefi}

类似的有频谱、相位角和功率谱如下：

\begin{align}
&|F(u,v)| =\sqrt{\left[ R^2(u,v)+I^2(u,v) \right]} \\
&\phi(u,v) = \arctan\frac{I(u,v)}{R(u,v)} \\
&P(u,v) = |F(u,v)|^2 =R^2(u,v)+I^2(u,v)
\end{align}


\subsection{程序代码说明}
因为有
\begin{equation}
\begin{split}
x(n)&=\frac{1}{N}\sum_{k=0}^{N-1}X(k)W^{-kn}_{N}=\frac{1}{N}\sum_{k=0}^{N-1}[X^*(k)W^{kn}_{N}]^*\\
&=\frac{1}{N}[\sum_{k=0}^{N-1}X^*(k)W^{kn}_{N}]^*\\
&=\frac{1}{N}(DFT[X^*(k)])^*
\end{split}
\end{equation}
即对已变换序列的共轭的DFT再共轭，可得原序列。所以在文件src/q4\_fft/fft.py中，最主要的部分就是1-D FFT正变换。1-D反变换，2-D正变换、反变换都是通过调用1-D FFT来实现的。

1-D FFT的计算步骤安照上节的16点FFT计算步骤图的方法进行。通过编写fft函数实现计算功能，在其中调用了辅助函数 \_fft\_one\_cal来进行上节图中一纵列的蝶形计算，不断对半减小参数unit\_size的值，直到为$2^1=2$完成计算。

ifft函数实现1-D 傅里叶反变换。fft2d实现2-D 傅里叶变换。ifft2d实现2-D 傅里叶反变换。

\subsection{实验结果与分析}
\subsubsection{算法性能}
既然是要求快速计算，所以将编写的FFT函数与不经优化的DFT函数和Numpy系统提供的高度优化的FFT函数进行比较。我们希望序列长度很大时，所编写的函数与系统提供的FFT函数最多不相差一个数量级。三种函数的比较如下：
\begin{longtable}{cccc}
\caption{算法性能} \\
%\toprule
%\multicolumn{3}{c}{算法性能} \\
%\midrule
%\endfirsthead
%\midrule
%\multicolumn{3}{c}{算法性能} \\
%\midrule
\endhead
\midrule
\multicolumn{3}{r}{接下页\dots} \\
\endfoot
\endlastfoot
\hline
样本序列长度&不经优化的DFT & 本程序FFT & Numpy的FFT \\
\hline
$2^8=256$&0.5387s&0.0116s&0.0002s\\
\hline
$2^{10}=1024$&8.4781s&0.0571s&0.0004s\\
\hline
$2^{16}=64K$&>1h&0.3922s&0.0315s\\
\bottomrule
\end{longtable}

可见，当序列长度不是太大时，简单的DFT与FFT差距不大，只差了一个数量级。估计这个时候的主要开销在Python环境中，算法的改进
影响很小，环境的改进就能有很大的收效。比如Numpy将FFT功能用C语言实现，其效果远好于用纯Python开发的版本。

当序列长度增长到$2^{16}=64K$时，简单的DFT函数已经慢的无法忍受了，但是使用了优化算法的FFT函数和Numpy的FFT函数的性能差距大大缩小了，只相差不到一个数量级了，并且几乎都是在瞬间完成的。

\subsubsection{图像结果}

\begin{longtable}{ll}
\caption{快速傅里叶变换} \\
%\toprule
%\multicolumn{2}{c}{快速傅里叶变换} \\
%\midrule
%\endfirsthead
%\midrule
%\multicolumn{2}{c}{快速傅里叶变换} \\
%\midrule
\endhead
\midrule
\multicolumn{2}{r}{接下页\dots} \\
\endfoot
\endlastfoot
\hline
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{tabular}{l}
频谱\\
\includegraphics[height=5cm]{../output/Q4_fft_amp.jpg}
\end{tabular}
&
\begin{tabular}{l}
相位谱\\
\includegraphics[height=5cm]{../output/Q4_fft_phanse.jpg} 
\end{tabular}
\\
\hline
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{tabular}{l}
反变换后的图像  \\
\includegraphics[height=5cm]{../output/Q4_ifft.jpg}
\end{tabular}
&
\begin{tabular}{l}
原图像与反变换的差图像 \\
\includegraphics[height=5cm]{../output/Q4_diff.jpg}
\end{tabular}
\\
\bottomrule
\end{longtable} 